Calculations of excess free energies of precipitates 
via direct thermodynamic integration across phase boundaries 
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We describe a technique for constraining macroscopic fluctuations in thermodynamic variables 
well-suited for Monte Carlo (MC) simulations of multiphase equilibria. In particular for multi- 
component systems this amounts to a statistical ensemble that implements constraints on both 
the average composition as well as its fluctuations. The variance-constrained semi-grandcanonical 
(VC-SGC) ensemble allows for MC simulations, in which single-phase systems can be reversibly 
switched into multiphase equilibria allowing the calculation of excess free energies of precipitates of 
complex shapes by thermodynamic integration. The basic features as well as the scaling and conver- 
gence properties of this technique are demonstrated by application to an Ising model. Finally, the 
VC-SGC MC simulation technique is used to calculate a /a' interface free energies in Fe-Cr alloys 
as a function of orientation and temperature taking into account configurational, vibrational, and 
structural degrees of freedom. 
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I. INTRODUCTION 

Multiphase systems are ubiquitous in nature. Since 
their properties are crucially dependent on the interfaces 
between the phases involved interface free energies con- 
stitute crucial input to kinetic and phase field models^ 
Their calculation, however, often constitutes a compli- 
cated taski^r— Here, an extension of conventional ther- 
mostatistical ensembles. Is described that enables calcu- 
lation of excess free energies via thermodynamic integra- 
tion across phase boundaries. The approach presented in 
this paper can be viewed as a generalization of the ex- 
tended Gaussian ensemble technique to multicomponent 
systems.— £ 

An equilibrium ensemble in statistical mechanics is a 
collection of macroscopic systems that have been pre- 
pared in the same thermodynamic state. This can be 
accomplished by constraining one of each pair of con- 
jugate thermodynamic variables that represent gener- 
alized displacements and forces, e.g. volume-pressure, 
entropy-temperature, or particle number-chemical po- 
tential. In single-phase regions of the phase diagram, 
there is a one-to-one correspondence between general- 
ized displacements and forces. This allows for straightfor- 
ward computation of free energy differences by thermo- 
dynamic integration. In contrast, in two-phase regions of 
the phase diagram, the mapping from generalized forces 
to generalized displacements becomes multi-valued. For 
instance the melting temperature maps to two distinct 
entropies, namely those of the solid S s and the liquid 
Si, respectively. Multiphase equilibria can be obtained 
by constraining the entropy to values S s < S < Si- If 
the conjugate forces are observables of the ensemble this 
enables thermodynamic integration across phase bound- 
aries. Conveniently, two important conjugate forces, 
namely temperature and pressure, are readily available 
in the microcanonical (NVE) ensemble. In contrast the 



chemical potential is not an observable of ensembles that 
control the number of particles or chemical composition 
in multicomponent systems. 

In this paper a technique is described for construct- 
ing ensembles, which constrain fluctuations of generalized 
displacements, that can be viewed as an extension and 
generalization of the Gaussian ensemble techniquei 7 -^ 
Such variance- constrained (VC) ensembles provide di- 
rect access to multiphase regions of the phase diagram 
while allowing for observation of both generalized forces 
and displacements. Thereby they enable thermody- 
namic integration across phase boundaries and compu- 
tation of free energies for compositions inside miscibil- 
ity gaps. In this paper this technique is employed to 
study precipitation in an Ising model system and a bi- 
nary Fe-Cr alloy. Through Monte Carlo simulations in 
the VC semi-grandcanonical (SGC) ensemble chemical 
potential-concentration isotherms in the shape of the 
van-der-Waals loop are obtained. In this fashion, one 
can extract interface free energies as a function of tem- 
perature and orientation. It is also demonstrated that 
interface free energies can be calculated by monitoring 
the growth and/or shrinkage of compact precipitates. 

The paper is organized as follows. After recapitu- 
lating key features of the related canonical and semi- 
grandcanonical ensembles, the next section introduces 
the VC-SGC ensemble and discusses strategies for sam- 
pling this ensemble using Monte Carlo (MC) simula- 
tions. The computation of interface free energies is 
demonstrated using a simple first-neighbor Ising model 
in Sect. MI A[ followed by an analysis of convergence and 
scaling properties in Sect. MI Bl Finally, in Sect. HVl the 
VC-SGC MC approach is applied to extract interface free 
energies in Fe-Cr alloys as a function of temperature 
and orientation while taking into account both config- 
urational and vibrational degrees of freedom. The paper 
is concluded in Sect. [Vj 
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II. ENSEMBLES FOR MULTI-COMPONENT 
SYSTEMS 

We start this section with a review of the canonical 
and semi-grandcanonical ensembles before deriving the 
variance-constrained semi-grandcanonical (VC-SGC) en- 
semble. We conclude by discussing practical aspects per- 
taining to sampling the VC-SGC ensemble using Monte 
Carlo simulations. 



A. Canonical and semi-grandcanonical ensembles 

Let us start by a brief discussion of multiphase equi- 
libria in multicomponcnt systems. For the sake of sim- 
plicity and clarity, we restrict ourselves to an immiscible 
binary alloy. The generalization to systems with an ar- 
bitrary number of species is straightforward. Consider 
a system of N particles confined in a box of volume V 
where each particle carries a spin of value or 1. An 
arbitrary configuration is specified by (x,<r), where x is 
a 3iV-dimensional vector describing the position of every 
particle in the box, and cr is an A-dimcnsional spin vec- 
tor. The number of spin 1 particles is n = $3 j=1 °i) an d 
their concentration is c = n/N. We denote the energy 
of a configuration by U{x,cr). The canonical partition 
function for this system at temperature T reads 



Z G (c,Af) 



n\(N -n)\ 



exp 



dx, (1) 



where = l/k B T and Af = {N, V, T} is the set of inde- 
pendent thermodynamic variables. Using the identity 
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we can rewrite the partition function in terms of an effec- 
tive potential where the compositional degrees of freedom 
have been integrated out: 

Zc (c, AT) = j-j f exp [-0F C (x- c, AT)) dx. (3) 
F c (as; c, AT) = -k B T\n S ( E * ~ n ) 

<T1—<TN \»=1 / 

x exp [-/3U"(a;,o-)] . (4) 

Fq (x; c, Af) is the effective potential energy landscape for 
the structural degrees of freedom of the multicomponcnt 
system at temperature T. In this way, we separate the 
compositional and topological degrees of freedom. In the 
following, we focus our discussion on the contribution of 
the compositional degrees of freedom to the free energy. 

We thus consider a frozen lattice of particles residing 
at x. For brevity we drop all references to x and denote 
the potential energy of a particular spin configuration as 



U(<r). The SGC ensemble represents a set of configura- 
tions that sample the spin degrees of freedom according 
to the Boltzmann distribution while also allowing the to- 
tal concentration to vary. To tune the average concen- 
tration an external chemical potential is applied to the 
system, which corresponds to modifying the potential en- 
ergy function as follows 



U s (cr; Afi) = U(a) + A^Ac(er) 



(5) 
(6) 



where A/i is the relative chemical potential that controls 
the average concentration and the free energy is given by 

F s (A/a, AT) = -fesTln ^ exp \-/3U s (a; An)] . (7) 



Inserting Eq. (JS|) into the above equation, the SGC par- 
tition function defined as Z§ = exp(— fg/fcsT) can be 
expressed in terms of the canonical free energy as 



Z s (An,Af) 



exp [-/? (F c (c, Af) + A/iNc)] dc. (8) 



The integrand in the above equation can be used to de- 
fine a concentration distribution function that is peaked 
around the average concentration (c) g . The condition of 
zero derivative at (c) s yields the well-known thermody- 
namic relation 



A, = -^«a> s ,A0 



(9) 



B. 



The variance constrained semi-grandcanonical 
ensemble 



Since for immiscible systems, one value of An maps 
to several compositions inside the miscibility gap stable 
multiphase coexistence cannot be established in the SGC 
ensemble. To overcome this restriction we modify the 
SGC ensemble in such a way as to control concentration 
fluctuations inside the miscibility gap. This is most eas- 
ily accomplished by adding a constraint that fixes the 
ensemble- averaged squared concentration (c 2 ). We call 
this the variance-constrained semi-grand canonical (VC- 
SGC) ensemble. 

This approach is analogous to umbrella sampling when 
calculating free energy barriers for reactions in molecular 
systemsi^^ In this method, given a particular reaction 
coordinate £, a harmonic external potential u b = 

K/2 (£ — £) is applied to bias the system toward posi- 
tions £ while at the same time the fluctuations arc re- 
strained by the force constant K. In the same spirit we 
can write the potential energy function of the VC-SGC 
ensemble as 



Uv(<t; <p, k) = U(a) + k (Nc(<t) + </>/2re) 2 



(10) 
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where and k are now thermodynamic variables that 
control the average concentration as well as its fluctua- 
tions. The harmonic external potential is parametrized 
such that Uy — > Us except for a constant shift when- 
ever k — y and <f> — > Afi. k can be given a physical 
intrepretation as the generalized force that controls con- 
centration fluctuations. This can be realized by putting 
the system in contact with a finite reservoir*^ Hence 
the VC-SGC ensemble, like the closely related extended 
Gaussian ensemble, is only applicable to finite systems. 
As all fluctuations vanish in the thermodynamic limit, so 
does k. It should thus diminish as the system size grows. 
This poses no difficulty in our study of interfacial free en- 
ergies since they also tend to zero in the thermodynamic 
limit. 

Analogously to the SGC ensemble [see Eqs. ([7]) and @] 
we can express the VC-SGC partition function in terms 
of the canonical free energy as 

Z w (<j>,K,M) = [ exp[-p(F c (c,N) + AU v )]dc, (11) 
Jo 

AUy = k (Nc + 4>/2k) 2 , (12) 

where A?7y corresponds to the harmonic bias potential 
of the umbrella sampling approach. The integrand in 
Eq. (fTTj) describes the probability distribution of the 
global concentration in the VC-SGC ensemble. It is a 
peaked function around the average concentration (c) v . 
The condition of zero derivative at (c) v yields the fol- 
lowing relation between (f> and k and the canonical free 
energy for the VC-SGC ensemble 

<f> + 2NK(c) v = -iV- 1 5Fc/9c((c) v ,A/'). (13) 

It is important to note that the right-hand side of 
Eq. (|13[) is the concentration derivative of the canoni- 
cal free energy at (c) v , which is a simulation observable. 
Hence no more unbiasing is necessary when employing 
the above technique for free energy integration. How- 
ever, for very small system sizes where the condition of 
zero derivative at (c) v may not hold, the simple proce- 
dure outlined above must be modified. The standard 
approach in such cases is to extract the probability dis- 
tribution of the system as a function of concentration via 
histogram methods^— and sampling bias must be cor- 
rected for as in the usual implementations of umbrella 
sampling ^ 10 ' 14 . 

Notwithstanding, the key advantage of the VC-SGC 
ensemble is that while the chemical driving force (right- 
hand side of Eq. (fT3| ) is an ensemble observable, the 
constraint on concentration fluctuations simultaneously 
allows multiphase equilibria to be stabilized. This is sim- 
pler and more efficient than the state-of-the-art methods 
for calculating concentration-derivatives of the free en- 
ergy within the canonical ensemble, which measure the 
work required to perform a gradual/instantaneous trans- 
mutation of one species into the other while keeping the 
particle number fix. While quite clever path-sampling 
algorithms with optimized estimators for calculating this 



work have been devised^ for the particular application 
discussed in this paper, the variance-constrained ensem- 
ble presents a simpler and more efficient route, since the 
free energy derivative can be obtained as an observable of 
the equilibrium ensemble, and no calculation of external 
work is required. Furthermore as discussed in detail in a 
previous publication^ the VC-SGC MC method allows 
for faster convergence to equilibrium in multiphase re- 
gions of the phase diagram than the canonical ensemble. 
Finally the VC-SGC ensemble is also very well suited for 
efficient parallel MC algorithms^ that enable simulations 
of very large systems containing millions of particles. 



C. Efficient sampling of the VC-SGC ensemble 

When carrying out MC simulations in the SGC en- 
semble according to Eq. ©, the parameter Afi is di- 
rectly proportional to the free energy derivative dF/dc, 
which renders choosing suitable parameters straightfor- 
ward. We will now demonstrate that choosing the rele- 
vant range of values for VC-SGC simulations is just as 
simple as in the case of the SGC ensemble. 

In the previous section the VC-SGC ensemble was de- 
rived in terms of the parameters 4> and k, which led to 
expressions that very much resembled the SGC ensemble. 
In practice it turns out that it is more convenient to use 
the substitutions k — R/N and <fi — R<fi. The VC-SGC 
potential defined above in Eq. (|10p then reads 

U v ((t;0,k) = U(tr) + RN(c + 0/2) (14) 

and the associated expression for the first derivative of 
the free energy [compare Eq. ([TB"]) ] is 

R (4> + 2 (c) v ) = -N~ l dF c /dc ((c) v , AO • (15) 

The above transformations effectively decouple the av- 
erage condtraint from the variance constraint parame- 
ter, which simplifies parameter selection for VC-SGC MC 
simulations as will be seen in the following. 

It is instructive to first consider a symmetric free en- 
ergy profile, in which case dF/dc = at the solubility 
limits where F is minimal as well as around 50% where 
F is maximal. From Eq. (|T5|) it follows that for non-zero 
R one can thus install a concentration of 50% by choosing 
<p = — 1. Furthermore one obtains (c) v — > 1 for <p w — 2 
and (c) v — s- for (f> m 0. These observations suggest a 
simple protocal for an efficient sampling of the full con- 
centration range: First choose a constant value for R and 
then vary <fi in the range —2.2 < cf> < 0.2. In our ex- 
perience this approach is very transferable and we have 
used it successfully in exactly this fashion for various sys- 
tems including for example the (very asymmetric) Fe-Cr 
system as described in Sect. ITVl 

In the following section it will be shown among over 
things that the results for the free energy derivative are 
insensitive to the variance constraint parameter R over a 
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FIG. 1. (Color online) (a) Derivative of free energy with re- 
spect to concentration at T = 1.5 ks for system size 8 x 8x 12. 
The gray areas indicate composition regions, in which spher- 
ical precitates and flat interfaces are formed. The bold green 
lines indicate fits to Eq. (|18[) . The inset shows the phase 
diagram for the Ising model described in the text, where tem- 
perature is measured in units of ks- (b) Free energy profile 
obtained by integration of data in (a). 



wide range of values. This provides for a lot of freedom 
in choosing hatn. 

It is straightforward to formulate a Monte Carlo 
algorithm^ for sampling the VC-SGC ensemble, where 
transmutation trial moves comprise 

1. selecting a particle at random 

2. flipping its spin (type) and 

3. computing the energy change AE and the concen- 
tration change Ac 

These trial moves are then accepted with probability 

Ay = min {l,exp [-f3(AE + kNc((j) + Ac + 2c))]} , 

(16) 

which satisfies detailed balance. 



III. BASIC FEATURES OF THE VC-SGC MC 
TECHNIQUE 

A. Interface free energies from direct 
thermodynamic integration 

In this section we illustrate the utility of the VC-SGC 
ensemble by studying phase segregation in a simple Ising 
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FIG. 2. (Color online) Interface free energies as a function of 
temperature and orientation for the Ising model described in 
the text. The inset shows a snapshot of a {100} interface from 
a VC-SGC MC simulation close to 50% composition (picture 
rendered with OVITO, Ref. [l7t) . 



model. We will calculate excess free energies of interfaces 
between two coexisting phases by direct thermodynamic 
integration from the single-phase region of the phase dia- 
gram into the miscibility gap. The Ising model is denned 
on a body-centered cubic (BCC) lattice with its Hamil- 
tonian parametrized as follows 



H 



9 SiSj, 



(17) 



where the summation runs over first-nearest neighbor 
pairs and Si G {—1,1}. The phase diagram of this system 
is shown in the inset of Fig. HJa). 

Our first goal in this section is to show that the VC- 
SGC ensemble, unlike the SGC ensemble, allows for cal- 
culation of free energy derivatives in both single and mul- 
tiphase regions of the phase diagram. To this end, sim- 
ulations were carried out using cells with basis vectors 
oriented along [100], [010], and [001]. In the VC-SGC 
simulations the constraint 4> (controlling average concen- 
tration) was varied from —2.05 to 0.05 in steps of 0.01 
at a constant variance constraint of k = 100. In case 
of SGC simulations the chemical potential difference A/x 
was changed between —3 and +3. At each value of <p 
(or A/i) the system was equilibrated for 5 x 10 3 MC 
sweeps, after which statistics were gathered for 15 x 10 3 
MC sweeps. 

The thus obtained relations between free energy 
derivative dF/dc and concentration c are shown in 
Fig. Ufa) f° r a temperature of 1.5 ks- The gap in the 
blue line in Fig. [TJa) corresponds to the miscibility gap. 
This reflects the inability of the SGC ensemble in sta- 
bilizing multiphase equilibria. The two concentrations 
delineating the boundary of the miscibility gap, i.e. 1% 
and 99% in Fig. [TJ are commonly referred to as binodals. 

Using the VC-SGC MC method F'{c) can be calcu- 
lated for all concentrations as shown by the yellow line 
in Fig. [TJa). Three regimes can be distinguished: At 
concentrations between the pure phases and the binodals 
(regime I), single-phase equilibria are obtained and the 
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FIG. 3. (Color online) Scaling of maximum of excess free 
energy [compare Fig. [IJb)] with inverse linear dimension for 
the Ising model described in the text. Simulation cell vectors 
are oriented along (100) directions inducing the formation of 
{100} interfaces. 



SGC-MC results are reproduced. In this regime, F'(c) 
has a positive slope. This behavior is preserved for com- 
positions between binodal and critical nucleation (ex- 
trema in dF/dc in Fig. []Ja)), which only support for- 
mation of subcritical minority-phase clusters (regime II) . 
In this regime there is a finite thermodynamic driving 
force for phase segregation but due to small sizes of the 
minority-phase clusters the intcrfacial free energy dom- 
inates and leads to their overall formation free energies 
being positive. Finally at concentrations beyond critical 
nucleation (regime III), supercritical clusters are formed 
and F'(c) is a monotonically decreasing function. In this 
regime the excess free energy increases due to the in- 
crease in the interface area of the cluster. Note that the 
van-dcr-Waals loop obtained in Fig. [TJa) is the equilib- 
rium solution under the constraint of finite system size. 
In particular, the concentration range spanned by regime 
II diminishes as the simulation cell increases, and even- 
tually vanishes in the thermodynamic limit. 

Studing the shapes of the supercritical clusters in 
regime III, we find several subregimes corresponding to 
different types of precipitates. For concentrations imme- 
diately beyond critical nucleation, minority-phase clus- 
ters are compact. The excess free energy (proportional 
to the interface area) as a function of concentration can 
thus be written as 



Aficfc) = £ [V co n(c - c )] 2/3 7 P 



(18) 



Here, £ is a shape factor, which equals (367r)3 for spheri- 
cal precipitates, and V ce n denotes the volume of the simu- 
lation cell. The bold dashed green lines in Fig. [T] indicate 
fits to Eq. (fT5|) and its concentration derivative. The ef- 
fective interface free energies calculated in this way are 
shown function of temperature in Fig. [2j 

At still higher concentrations, the precipitates become 
so large that spherical clusters can no longer be contained 
inside the finite-sized simulation cell and a transition to 
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FIG. 4. (Color online) Standard deviation of concentration 
at a temperature of T = 2 fcsas a function of (a) variance con- 
straint parameter and (b) system size. Data in (b) measured 
at 50% composition. The black dashed line indicates scaling 
with \j\fx. 



a cylindrical shape induced by periodic boundary condi- 
tions is observed. While the free energy profile recorded 
in this region is of limited physical interest per se, it pro- 
vides the continuous derivative needed for integration of 
the full excess free energy. 

Finally in the concentration range around 50% in 
Fig.QJa), the excess free energy assumes a constant value 
due to formation of supcrlatticcs consisting of alternating 
slabs of the two phases, which follow the periodicity of 
the simulation cell. The excess free energy at these con- 
centrations corresponds to the free energy cost associated 
with two flat interfaces with the free energy density given 
by 



7fla t = AF XC /A, 



(19) 



where A is the cross-sectional area of the computational 
cell. By changing the geometry of the simulation cell it 
is possible to change A and thereby interface orientation. 
To see this remember that during the course of a simula- 
tion the system will establish the interface configuration 
that minimizes AF XC not 7fl a t- Since the cross section A 
for different orientations can be changed by varying shape 
and size of the simulation cell, it is possible to stabilize 
different interface orientations as exemplified in Fig. [2] 
To obtain for example the free energy density of {111} 
interfaces, we used an orthorhombic simulation cell with 
lattice vectors oriented along [110], [112], and [111] and 
2x3x8 unit cells. In this configuration the excess free 
energy AF C for {100} and {110} interfaces is larger than 
for {111} interfaces. Similarly shortening the dimensions 
perpendicular to [110] and extending the cell parallel to 
this direction stabilizes a {110} interface. In this fashion 
it is possible to not only extract the temperature but also 
the orientation dependence of the interface free energy, 
which leads to the data displayed in Fig. [5] 
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To verify the reliability of the thus computed interfa- o 
cial free energies using the VC-SGC MC technique, we 

have also studied the dependence of the excess free en- " 1 

ergies of the two-phase equilibria inside the miscibility .> 

gap on the size of the simulation cell. Since the excess > 

free energy is dominated by the interfacial free energy it -g 

should scale inversely with the longest linear dimension g; 

of the simulation cell. This behavior is indeed observed £j _., 

as shown in Fig. [3j % 
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B. Convergence and parameter sensitivity 

In this section we study the convergence properties of 
VC-SGC MC simulations and the sensitivity of the cal- 
culated free energies to the choice of the variance con- 
straint parameter. Figure 0] shows the dependence of the 
standard deviation of concentration on (a) variance con- 
straint parameter and (b) system size. In both cases the 
standard deviation scales inversely with the square root 
of the respective parameter. To understand this behav- 
ior consider again Eq. (fTTj). The standard deviation of 
concentration in the VC-SGC ensemble can be expressed 
as the second moment of its probability distribution as 
follows 

1 r 1 

s 2 v = — / c 2 exp [— jS (F c (c, TV) + At/*)] dc - (c) v . 

Jo 

(20) 

For large enough system sizes, the ensemble prob- 
ability distribution becomes a normal distribution 
^A/nexp (— A(c — -B) 2 ), with its mean B = (c) v , and 
its width A inversely proportional to the standard devi- 
ation: Sy = 1/2 A. Now the Gaussian width A is related 
to the second derivative of the distribution taken at the 
average concentration. We thus obtain the following re- 

? 2 
3 v 



1 



d*F c ({c) v ,N) 
dc 2 



2RN 



(21) 



This is a very important relation that clearly predicts 
the dependence of the standard deviation of concentra- 
tion in the VC-SGC ensemble on the variance constraint 
parameter R as well as on the system size N as observed 
in Figs.SJa) andQJb). It also provides a clear picture of 
why and under which conditions the VC-SGC ensemble 
can stabilize multiphase equilibria. Equation (|21j) illus- 
trates that whenever R — > as in the SGC ensemble, sy 
diverges for values of (c)y, for which the second deriva- 
tive of the free energy function Fc (c, A/") becomes nega- 
tive. For these compositions, stable single-phase equilib- 
ria cannot be achieved. Introducing a non-zero variance 
constraint R leads to finite sy. However, in order for R 
to stabilize compositions inside the miscibility, it has to 
be chosen large enough to make the right-hand side of 
Eq. (f2"Tj) positive. Hence for practical calculations, the 
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FIG. 5. (Color online) Effect of variance constraint parame- 
ter R on (a) free energy derivative and (b) standard deviation 
of concentration at T — 2 ks for system size 4x4x8. 



best choice for R is a value slightly larger than the max- 
imum of l^c ( C ' J ^0I ms ide the miscibility gap. Ac- 
cording to Fig. 01a) this threshold is about R = 20 for a 
system with 4x4x8 unit cells and close to 10 for a size 
of 8 x 8 x 8 unit cells. 

For k — y the VC-SGC ensemble approaches the SGC 
ensemble, for which the standard deviation should di- 
verge for concentrations inside the miscibility gap, a be- 
havior that is in fact clearly visible in Fig. 2J The figure 
also explicitly demonstrates that in the limits R — > oo 
and/or N — > oo concentration fluctuations go to zero 
and thus approach the canonical ensemble. 

Figure [5^a) illustrates the important fact that the free 
energy derivative obtained via Eq. (|13[) is indeed unaf- 
fected by the strength of the variance constraint over a 
wide range of parameter values. The lower limit of the 
acceptable parameter range is set by the divergence of 
the standard deviation that was described in the previ- 
ous paragraph and that is shown as a function of concen- 
tration in Fig.(5Jb). The upper limit, however, is due to 
decreasing acceptance probability of the MC trial moves, 
which as illustrated in Fig.[S]diminishcs only slightly with 
R over several orders of magnitude before dropping more 
rapidly. As a result for large values of R one can no longer 
gather sufficient statistics to compute meanigful thermo- 
dynamic averages. The range within which R can be se- 
lected nonetheless extends at least over two to three or- 
ders of magnitude depending on system size. Since large 
R values lead to a reduction of the acceptance probabil- 
ity, it is, however, recommended to choose R as small as 
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FIG. 6. (Color online) Effect of variance constraint parame- 
ter R on (a) free energy derivative and (b) standard deviation 
of concentration at T = 2 fcs for system size 4x4x8. 



possible while maintaining an acceptable standard devi- 
ation. Figure @] suggests that a reasonable value is about 
1%. 



IV. INTERFACE FREE ENERGIES IN THE 
FE-CR ALLOY SYSTEM 

Up until now we have considered a simple Ising model. 
It can be considered as the most simple representative 
of alloy cluster expansions ; 18 ' 19 which are widely used 
to study complex multi-component alloys. The VC-SGC 
MC technique is, however, equally applicable to empirical 
potential models and then enables computation of excess 
free energies taking into account not only configurational 
but also structural and vibrational degress of freedom^ 
To illustrate this point, we now consider the calculation 
of excess free energies in BCC the Fe-Cr binary alloy sys- 
tem as described by the concentration-dependent embed- 
ded atom methodi 21 ' 22 The latter reproduces the phase 
diagram below the Curie temperature featuring both sub- 
stantial solubility of Cr in Fe (cv-phase) and a pronounced 
miscibility gap as shown in the inset of Fig. [7ta). 

Simulations were carried out using orthorhombic cells 
containing between 384 and 480 atoms depending on 
sample orientation. The variance constraint parameter 
was set to K=100cV/atom and <j> was varied between 
—2.1 and 0.1 (see Sect. Ill Cj) . Structural relaxations and 
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FIG. 7. (a) Derivative of free energy with respect to con- 
centration obtained from MC simulations in the SGC and 
VC-SGC ensembles at 1000 K. The inset shows the calcu- 
lated phase diagram, which is in excellent agreement with the 
published phase diagram of the potential used in the present 
work.— (b) Free energy profile obtained by integration of the 
data in (a) and (c) the corresponding excess free energy. The 
central dark gray area indicates the region within which flat 
interfaces are obtained while the dark gray areas to the left 
and right show the regions in which compact precipitates are 
obtained. 



thermal vibrations were sampled using displacement MC 
trial moves, while volume trial moves were employed to 
sample zero pressure conditions. The system was evolved 
for 6,000 MC steps at each value of (f>. Statistics were 
gathered over the last 5,000 steps. 

First we mapped out dF/dc using SGC MC simula- 
tions as shown for a temperature of 1000 K by the bold 
blue line in Fig. [TJa). The gap in the blue line corre- 
sponds to the miscibility gap, which again reflects the in- 
ability of the SGC ensemble to stabilize multiphase equi- 
libria. Using the VC-SGC MC method we then calcu- 
lated OF/ dc over the entire concentration range as shown 
by the yellow line in Fig. EJa). The free energy profile 
obtained in this fashion again coincides with the SGC 
results in the single phase regions and inside the misci- 
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FIG. 8. Interface free energies for Fe-Cr model alloy for dif- 
ferent interface orientations as well as spherical precipitates. 



bility gap shows the same features that were discussed in 
Sect. MI XI for the Ising model. 

Once again it is possible to extract free energy densi- 
ties corresponding to spherical and flat interfaces from 
the excess free energy curves, as shown in Fig. [Tfc) . The 
results of this analysis are summarized in Fig. [8] Note 
that in this case the computation of zero-K interface en- 
ergies is significantly more complicated than in the case 
of the Ising model. The asymmetric phase diagram fea- 
tures a large solubility of Cr in a-Fe that remains finite as 
temperature goes to zero and is accompanied by short- 
range ordering.— We have therefore not attempted to 
construct a zero-K a /a' interface. 

Two features in Fig. [S] are particularly noteworthy in 
comparison with Fig. [5] The ordering of the interface 
orientations is different with {100} interfaces now being 
lowest in energy, and the anisotropy is much smaller. 

It is furthermore instructive to separate the different 
contributions to the interface free energy. First, by re- 
peating the calculations for {100} interfaces without dis- 
placement trials moves^ it is possible to separate vibra- 
tional and configurational degrees of freedom. The result 
of this analysis is shown in Fig.[9ja), which illustrates the 
magnitude of the vibrational contribution to the interface 
free energy. Since the simulations also provide average 
internal energies one can furthermore separate explicitly 
the interface internal energies and entropies. As shown 
in Fig. |HT,b) and (c) the internal energies are hardly af- 
fected if vibrations and atomic relaxations (other than 
volume changes) are suppressed. Instead their contribu- 
tions show up as a practically temperature independent 
contribution to the entropy. Figures IHta) and (b) finally 
show that interface free and internal energies extrapolate 
to a value of about 372 ± 2mJ/m 2 at zero K. 

V. CONCLUSIONS 

In this paper we have introduced a free energy integra- 
tion technique along compositional degrees of freedom 
that allows for efficient calculations of free energies of 
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FIG. 9. (a) Free energy, (b) internal energy, and (c) entropy 
of {100} interface in Fe-Cr model alloy. The dashed lines 
represent fits to second degree polynomials. 



multiphase systems with various degrees of heterogene- 
ity. The construction of the underlying VC-SGC ensem- 
ble represents a generalization of the extended Gaussian 
ensemble technique to multicomponent systems 

To illustrate the basic features of the VC-SGC MC 
technique as well as aspects such as scaling and con- 
vergence, an extensive characterization of a first-nearest 
neighbor Ising model was carried out. In this way it was 
shown that the VC-SGC MC technique allows (z) sta- 
bilizing arbitrary concentrations inside (and outside) the 
miscibility gap and (ii) simultaneously computing deriva- 
tives of the free energy with respect to concentration. 
Efficient sampling of the entire concentration range is 
straightforward since in essence it merely requires setting 
the variance parameter R, which can be chosen rather de- 
liberately. 

Since VC-SGC MC simulations provide free energy 
derivatives as a continuous function of composition, ther- 
modynamic integration can be carried out along the con- 
centration axis, which allows determination of interface 
free energy densities. Specifically by monitoring the free 
energy derivative as precipitate size and shape evolve, it 
is possible to extract the excess free energies of compact 
precipitates. By pushing the system close to 50-50 phase 
composition while varying simulation cell shape and size, 
one can also derive the interface free energies of flat in- 
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tcrfaccs as a function of orientation. 

To demonstrate the applicability of the VC-SGC ap- 
proach in more general cases, we finally considered Fe-Cr 
alloys as described by an empirical potential scheme. The 
interface free energies for this system were determined 
taking into account not only the configurational but also 
vibrational and structural degrees of freedom. Our re- 
sults indicate a strong temperature but weak orientation 
dependence of the interface free energies, suggesting that 
both a and a 1 precipitates should adopt a spherical shape 
at practically all relevant temperatures. 

As described in Ref. M, the VC-SGC MC ensemble 
leans itself to efficient parallelization, enabling simula- 
tions of precipitation in multi-million atom samples^ It 
should finally be stressed that our approach is also di- 
rectly applicable for obtaining interface free energies and 



nucleation barriers from cluster-expansions ; 18 ' 19 which 
are widely applied to study complex multi-component 
alloys. 
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